!
! Copyright (C) 2000-2013 C. Hogan and the YAMBO team 
!              https://code.google.com/p/rocinante.org
! 
! This file is distributed under the terms of the GNU 
! General Public License. You can redistribute it and/or 
! modify it under the terms of the GNU General Public 
! License as published by the Free Software Foundation; 
! either version 2, or (at your option) any later version.
!
! This program is distributed in the hope that it will 
! be useful, but WITHOUT ANY WARRANTY; without even the 
! implied warranty of MERCHANTABILITY or FITNESS FOR A 
! PARTICULAR PURPOSE.  See the GNU General Public License 
! for more details.
!
! You should have received a copy of the GNU General Public 
! License along with this program; if not, write to the Free 
! Software Foundation, Inc., 59 Temple Place - Suite 330,Boston, 
! MA 02111-1307, USA or visit http://www.gnu.org/copyleft/gpl.txt.
!
subroutine surfdriver(Xen,Xk,q,wv,X) ! remove q
    
  use pars,                    only : SP, PI, schlen, lchlen
  use units,                   only : HA2EV
  use com,                     only : msg
  use memory_m,                only : mem_est
  use parser_m,                only : parser
  use optcut,                  only : loptcut, setup_cutoff, print_cutoff, dc
  use electrons,               only : levels
  use R_lattice,               only : bz_samp, bare_qpg,q_norm
  use IO_m,                    only : io_control, OP_RD_CL, OP_WR_CL, VERIFY, REP, DUMP
  use X_m,                     only : X_t, X_epsilon, X_mat, X_alloc
  use frequency,               only : w_samp 
  use optcut,                  only : loptcut, DIP_q_dot_iR_cut, setup_cutmat ! just for dealloc
  use vec_operate,             only : v_norm
  use parallel_m,              only : pp_redux_wait
  use wave_func,               only : WF_free
  use ras_module
  use surface_geometry
  use bulkeps,                 only : lbulkerr, GetBulkEps
  use eels_kinematics
  use eels_detector
  use model_loss_function,     only : lossformlg, lossfm, print_eels_form
  implicit none
  type(levels),       intent(in) :: Xen    ! Energies
  type(bz_samp),   intent(inout) :: Xk, q     ! K points
  type(X_t),       intent(inout) :: X      ! Polarization
  type(w_samp)                   :: wv     ! Energy ranges
!ws
  integer                        :: fr(2),iq, ipol, iw
  real(SP)                       :: qpard(2), q0v(3,3)
  character(schlen)              :: lch

! Full spectrum arrays 
  complex(SP), allocatable       :: eps_b(:), eps_i(:,:), surf_i(:,:)
  real(SP),  allocatable         :: lossf(:,:)

  integer                        :: nw, npol, io_err, ID
  logical                        :: lfail, lsurferr, lreelserr, lraserr
  character(lchlen)              :: bulkmsg

  call section('*','Surface spectroscopy module')

  lbulkerr = .false.; lsurferr = .false.; lraserr = .false.; lreelserr = .false.
  lfail = .false.
  !
  ! Polarizations
  !
  if(.not.lras.and..not.lreels) return 
  npol = 2
  if(lreels) npol = 3
  !
  ! Surface geometry and spectral analysis
  !
  call section('+','Surface geometry and spectral analysis')
  call setup_surface_geometry( lfail )
  call print_surface_geometry
  call setup_cutoff( lfail )
  call print_cutoff
  if(loptcut) then
    call setup_gvecaff
    call setup_cutmat
  endif
  call setup_loc(X%ib(1), Xen%nbm, Xen%nbf+1, X%ib(2), Xk%nibz)
  !
  ! REELS input parameters
  !
  if(lreels) then
    call section('=','REELS setup')
    call setup_eels_kin( lfail )
    if(lfail) lreelserr = .true.
    call print_eels_kin
    call setup_eels_det( lfail)
    if(lfail) lreelserr = .true.
    call print_eels_det
    call check_eels_det( lfail )
    if(lfail) lreelserr = .true.
    call print_eels_form
    call print_eels_geometry
    call print_eels_general( npol )
  endif
  q0v(:,1) = q0x ; q0v(:,2) = q0y ; q0v(:,3) = norm
  !
  ! Energy grids and energy dependent arrays
  !
  call FREQUENCIES_setup(wv) 
  call X_pre_setup(Xen,X)
  nw = wv%n(1) ! Have to use freqs_setup in case wv%n(1) is lowered (%)

  allocate( eps_i(nw, npol) )  
  allocate( eps_b(nw) )
  allocate( lossf(nw,4) )
  !
  ! Read and process the bulk data
  !
  call section('=','Bulk epsilon data import')
  call GetBulkEps(real(wv % p), nw, eps_b, bulkmsg)
  if(lbulkerr) lreelserr = .true.
  if(lbulkerr) lraserr = .true.
  call msg('nrs',bulkmsg)
  if(lbulkerr) call msg('nrs','Calculation will continue, but '//&
&   'quantities requiring the bulk epsilon will not be computed.')
  !
  ! Compute epsilon
  !
  call section('=','Calculating eps cell.')

  allocate( X_epsilon(2,nw) ) 
  !
  ! Main loop over polarizations
  !
  do ipol = 1, npol

    iq = 1

    select case(ipol)
    case(1)
       X%q0 = q0x
       qname = 'x'
    case(2)
       X%q0 = q0y
       qname = 'y'
       q0v(:,2) = q0y
    case(3)
       X%q0 = norm
       qname = 'z'
       q0v(:,3) = norm
    end select
    write (lch,'(3a,3f10.5)') &
&            'Absorption @ q || ',qname,' : ', X%q0/v_norm(X%q0)
    call msg('ns',trim(lch)) ! <-- to std out + newline
       
    call X_alloc('X',(/X%ng,X%ng,wv%n(2)/))

    X_epsilon = (0.,0.)

    do iw = 1, wv % n(1), wv % n(2) ! loop over blocks
       
      fr = (/ iw, iw+wv%n(2)-1 /) ! the small block of frequencies
       
      call X_os(X_mat,iq,fr,Xen,Xk,wv,X)
       
      if(loptcut) X_mat(1,1,:) = X_mat(1,1,:)*az/dc ! Scale if cutoff present

      X_epsilon(1,fr(1):fr(2))=1.0_SP-X_mat(1,1,:)*4.0_SP*pi/q_norm(iq)**2
       
    enddo 
    !
    ! Clean up for each q
    !
    call X_alloc('X') 
    !
    ! Copy to slab eps. array
    !
    do iw = 1, nw
      eps_i(iw,ipol) = X_epsilon(1,iw)
    enddo
    call X_alloc('OptOsc')  ! DIP_iq_dot_r
    if(loptcut) deallocate(DIP_q_dot_iR_cut)
     
  enddo    
  deallocate(X_epsilon)
  !
  ! Dump the eps_i (cell) to file for plotting
  !
  call section('=','Write eps cell for plotting')
  call write_eps_tensor( q0v, eps_i, real(wv % p), nw, npol, 'cell' ) 
  !
  ! Find slab symmetries
  !
  allocate( surf_i(nw, npol ) )
  call section('=','Slab symmetries ')
  call get_slab_symmetry( loptcut) 
  !
  ! Extract surface dielectric function from cutoff/bulk subtraction
  !
  call section('=','Surface dielectric function')
  call extract_eps_surf(surf_i, eps_i, npol, eps_b, nw, lbulkerr, lsurferr) 
  if(lsurferr) lraserr = .true. ;  if(lsurferr) lreelserr = .true.
  if(.not.lsurferr) call write_eps_tensor( q0v, surf_i, real(wv % p), nw, npol, 'surf' ) 
  !
  ! Calculate RAS and REELS in 3-layer method
  !
  if (lras.and..not.lraserr) then
    !
    call section('=','RAS')
    call CalcRAS(real(wv % p), surf_i(1:nw,1), surf_i(1:nw,2), &
&                              eps_b(1:nw), nw, d_surf,'-rpa')
    !
  endif
  !
  if(lreels.and..not.lreelserr) then
    
    lch = 'REELS in model form: '//trim(lossformlg(lossfm))
    call section('=',trim(lch) )
    call CalcREELS(lossf, real(wv % p), surf_i, eps_b, nw, qpard)
    call write_reelsdiff(lossf, real(wv % p), nw ) 
    !
  endif
  !
  ! Clean up
  !
  deallocate( eps_i, eps_b, surf_i )
  deallocate(wv%p)
  call mem_est("W-p")
  deallocate(lossf)
  deallocate(Xk%pt)
  call mem_est("k_pt")
  call WF_free()
  call ras_clean

  call pp_redux_wait

  return
end subroutine surfdriver

subroutine ras_clean
 use memory_m,      only:mem_est
 use bulkeps,        only : clean_bulkeps
 use surface_geometry, only : gvecaff
 use optcut, only : end_optcut,loptcut
  use frequency,     only:rg_index_bg,bg_npts,cg_pt,cg_index_bg
  use X_m,                     only : X_poles_tab
  use R_lattice,        only: g_vec, bare_qpg, g_rot, qindx_free
  use D_lattice,        only: dl_sop
  use wave_func,        only: wf_nc_k, wf_igk

  call qindx_free
  call clean_bulkeps
  if(allocated(gvecaff)) deallocate(gvecaff)
  if(loptcut) call end_optcut

  deallocate(g_vec,dl_sop,wf_nc_k,wf_igk)
  call mem_est("g_vec dl_sop wf_nc_k wf_igk")

  if(allocated(X_poles_tab)) then
    deallocate(X_poles_tab)
    call mem_est(" X_poles_tab ")
  endif

  if(allocated(rg_index_bg)) call mem_est("RGi BGn CGp CGi")
  if(allocated(rg_index_bg)) deallocate(rg_index_bg)
  if(allocated(bg_npts)) deallocate(bg_npts)
  if(allocated(cg_pt)) deallocate(cg_pt)
  if(allocated(cg_index_bg)) deallocate(cg_index_bg)

  deallocate(bare_qpg)
  call mem_est("bare_qpg")

  if(allocated(g_rot)) deallocate(g_rot)

  return
end subroutine ras_clean


subroutine ras2close
 use memory_m,      only:mem_est
  use R_lattice,        only: ng_in_shell,E_of_shell
  if(allocated(E_of_shell)) then
    deallocate(E_of_shell)
  endif
  if(allocated(ng_in_shell)) deallocate(ng_in_shell)

end subroutine ras2close
 
subroutine print_eels_general( npol )
  use com, only : msg
  implicit none
  integer, intent(in)      :: npol

  if(npol.eq.2) call msg('r',&
&   'Double q-polarization and relative difference calculations [RAS].')
  if(npol.eq.3) call msg('r',&
&   'Triple q-polarization calculation selected [RAS, HREELS].')
  return
end subroutine print_eels_general

subroutine setup_loc(iv1,iv2,ic1,ic2)
  use com, only : msg
  use pars,                    only : SP
  use IO_m,                    only : io_control, OP_RD_CL, OP_WR_CL, VERIFY, REP, DUMP
  use ras_module,              only : loctype, locfac, transloc, upperlim, lowerlim
  use electrons,               only : n_bands
  use R_lattice,               only : nkibz
  implicit none
  integer, external :: io_loc
  integer :: ID_loc, io_err, ik, iv, ic
  integer, intent (in) :: iv1, iv2, ic1, ic2
  real(SP)                        :: loc(n_bands, nkibz)

  allocate(transloc(ic2,iv2, nkibz))
  transloc = .true.
  if(loctype.eq."no") return
  call io_control(ACTION=OP_RD_CL, COM=REP, SEC=(/1/), ID=ID_loc)
  io_err = io_loc(loc,upperlim, lowerlim,ID_loc)
  if(io_err.ne.0) return

  call msg('rs','Restricting transitions to type: ',loctype)
  do ik = 1, nkibz
    do iv = iv1, iv2
      do ic = ic1, ic2
        if(loctype.eq."ss") transloc(ic,iv,ik) = &
&                           (loc(iv,ik) >= locfac).and.(loc(ic,ik) >= locfac)
        if(loctype.eq."sb") transloc(ic,iv,ik) = &
&                           (loc(iv,ik) >= locfac).and.(loc(ic,ik) <  locfac)
        if(loctype.eq."bs") transloc(ic,iv,ik) = &
&                           (loc(iv,ik) <  locfac).and.(loc(ic,ik) >= locfac)
        if(loctype.eq."bb") transloc(ic,iv,ik) = &
&                           (loc(iv,ik) <  locfac).and.(loc(ic,ik) <  locfac)
!DEBUG>
!write(88,'(3i4,l2,2f7.3,1x,f7.3)') ik,iv,ic,transloc(ic,iv,ik),loc(iv,ik),loc(ic,ik),locfac  ! DEBUG
!DEBUG<
      enddo
    enddo
  enddo
  return
end subroutine setup_loc
